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ABSTRACT 

The vertical shear instability (VSI) offers a potential hydrodynamic mechanism for angular momen¬ 
tum transport in protoplanetary disks (PPDs). The VSI is driven by a weak vertical gradient in the 
disk’s orbital motion, but must overcome vertical buoyancy, a strongly stabilizing influence in cold 
disks, where heating is dominated by external irradiation. Rapid radiative cooling reduces the effec¬ 
tive buoyancy and allows the VSI to operate. We quantify the cooling timescale tc needed for efficient 
VSI growth, through a linear analysis of the VSI with cooling in vertically global, radially local disk 
models. We And the VSI is most vigorous for rapid cooling with U < — 1) in terms of 

the Keplerian orbital frequency, Hk; the disk’s aspect-ratio, /i <C 1; the radial power-law temperature 
gradient, q; and the adiabatic index, 7 . For longer U, the VSI is much less effective because growth 
slows and shifts to smaller length scales, which are more prone to viscous or turbulent decay. We 
apply our results to PPD models where is determined by the opacity of dust grains. We And that 
the VSI is most effective at intermediate radii, from ^ 5AU to ^ 50AU with a characteristic growth 
time of ^ 30 local orbital periods. Growth is suppressed by long cooling times both in the opaque 
inner disk and the optically thin outer disk. Reducing the dust opacity by a factor of 10 increases 
cooling times enough to quench the VSI at all disk radii. Thus the formation of solid protoplanets, a 
sink for dust grains, can impede the VSI. 


1. INTRODUCTION 

Understanding how disks transport mass and an¬ 
gular momentum underlies many problems in astro¬ 
physi cs, including star and planet formation (lArmitagel 
l2011li . The turbulence associated with many transport 
mechanisms is particularly im portant for dust evolution 
and planetesimal form ation (lYoudin fc Lithwicld I2007t 
iChiang fc Youdinll20ir)ll . 

Magneto-hydrodynamic (MHD) turbulence 
driven by the magne to-rotational instability (MRI, 
iBalbus fc Hawle^IlQQlf) has long been the most promis¬ 
ing transport mechanism in low mass disks with weak 
self-gravity. However, many parts of protoplanetary 
disks (PPDs) are cold, have lo w levels of ionization , 
and do not support the MRI (|Blaes fc: BalbusI II994 
iSalmeron fc Wardig I2003D . Recent simulations suggest 
that significant portions of PPDs fail to develop MHD 
turbulence le.g.lSimon et al.ll2013l: iLesur et'n]l20l4 iBail 
[W^ICre.ssel et al. 11291,^ ! 

A purely hydrodynamic mechanism could circumvent 
difficulties with non-ideal MHD, but must overcome the 
strong centrifugal stability imposed by the positive ra¬ 
dial specific a ngular momentum gradient in nearly Ke¬ 
plerian disks (IBalbus et al.lll99'^ . One possible route 
to hydrodyna mic turbulence is the vertical shear in¬ 
stability Ybl. lUrpin fc Brandenbu7^ll998l : IUrpinl[2?)0^ 
INelson et al.ll201^ hereafter INI311 . The basic mechanism 
of the VSI in disks was firs t identified in the context of 
differenti ally ro t ating stars llGoldreich fc SchubertllI967L 
hereafter IGS671 lFrickelll968ll . The VSI arises when ver¬ 
tical shear, i.e. a variation in orbital motion along the 
rotation axis, destabilizes inertial-gravity waves, which 
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are oscillations with rotation and buoyancy as restoring 
forces. Vertical shear occurs wherever the disk is baro- 
clinic, i.e. when constant density and constant pressure 
surfaces are misaligned. Baroclinicity, and thus vertical 
shear, is practically unavoidable in astrophysical disks, 
except at special locations like the midplane. 

To overcome centrifugal stabilization, the VSI triggers 
motions which are vertically elongated and radially nar¬ 
row. Vertic al elongation taps the free energy of the ver¬ 
tical shear (|Umurhan et al.l[2013fl . but is also subject to 
the stabilizing effects of vertical buoyancy if the disk is 
stably stratified. To overcome vertica.1 buo yancy, the VSI 
requires a short cooling time (IGS67I : IN13I1 . Rapid radia¬ 
tive cooling, i.e. a short thermal relaxation timescale, 
brings a moving fluid element into thermal equilibrium 
with its surroundings, thereby diminishing buoyancy. An 
isothermal equation of state implies instantaneous ther¬ 
mal r elaxation, and is t h e ideal context for studyi ng the 
VSI (lUrpid 20031: INI3I: iMcNallv fc Pessahl 12014 here¬ 
after IMP 14 Barker fc Lattei]l2015l hereafter lBLlSll . 

Alternatively, vertical buoyancy can be eliminated by 
strong internal heating, i.e. by the onset of convection, 
so that the disk becomes ne utrally stratified in the ver¬ 
tical direction (IN13I : iBLlSfl . However, realistic PPDs 
should be vertically stably stratified in the outer re¬ 
gions, beyo nd ^ 1-5 AU, wh ere heating is irradiation 
dominated (jBitsch et al.ll2015l) . Even in the inner disk, 
strong vertical buoyancy is possible if accr etion heating is 
weak or is concentrate d in surface layers ()Gammielll996l : 
iLesniak fc Deschll2011l) . 

Understanding the VSI in real disks therefore requires 
considering finite, non-zero cooling times, U- Non-linear 
hydrodynamical simulations with a prescribed U find 
that VSI turbulence in stably stratified disks requires 
rapid cooling with tc shorter than orbital timescales 
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(|N13D . When the cooling time is short enough, the VSI 
can drive moderately strong transport, with Re ynold s 
stresses of a ~ 10“^ times the mean pressure (INIST ). 
Simulations with realistic radiative transfer, in lieu of 
a fixed tc, find that the VSI in irradiated disks drive 
transport with a ^ 10“'^ in a ~ 2 - lOAU PPD model 
(IStoll fc KlevlIMl . 

Studying the linear growth of the VSI is neces¬ 
sary for understanding how it can ultimately drive 
turbulent transport. The pioneering linear analy- 
ses of the VSI considered verticallv local disturbances 
(|IJrDin fc Brandenburd 119981 lIJroii] I2003I1 . A vertical 
global analysis (e.g. ' INLH IMP14 IBLl,»ill is essential for 
understanding how vertically elongated disturbances in¬ 
teract with the disk’s vertical structure. Moreover a ver¬ 
tically global analysis allows more direct comparisons 
with modern numerical simulations. This work gener¬ 
alizes previous vertically global analyses by including a 
finite cooling time. 

This paper is organized as follows. In ^we explain, 
without derivation, our main result for the critical cool¬ 
ing timescale, beyond which VSI growth is suppressed. 
We develop our disk model in ^ and derive linearized 
perturbation equations in m Section [5] contains our 
main analytic results, leading to the derivation of the 
critical cooling time in 115.31 In ^ we analyze linear 
VSI growth, and numerically confirm the critical cooling 
time. We apply our results to PPDs in ^ with cooling 
times derived from dust opacities. We discuss caveats 
and extensions in ^ We conclude in ^ with a sum¬ 
mary. Some technical and background developments are 
explained in the appendices. 


2. WHY MUST COOLING BE SO FAST? 

In a stably stratified thin disk, the VSI requires a ther¬ 
mal timescale tc significantly shorter than the disk dy¬ 
namical timescale, 

tc < ( 1 ) 


where Hk is the Keplerian frequency (N13). This require¬ 
ment — that the cooling time be much shorter than an 
orbital period, which in turn is much smaller than the 
relevant oscillation period of vertically elongated grav¬ 
ity waves — is quite stringent. It highlights the fact that 
vertical buoyancy is strongly stabilizing. Rapid radiative 
damping is thus required to weaken buoyancy and allow 
the weak vertical shear to drive instability. 

To roughly quantify the required smallness of tc, we 
consider a vertically isothermal disk with aspect-ratio h, 
radial temperature profile T oc r'J and adiabatic index 
7 > 1 so the disk is stably stratified. For a PPD, h ^ 
0.05, q ^ —0.5 and 7 ~ 1.4. 

For a thin disk {h <C 1), the vertical shear rate varies 
with height, z, from the midplane as 


dVt 


r 


dz 


^ (^) , ( 2 ) 


where H — hr is the characteristic disk scale height. 

This destabilizing shear competes with the stabilizing 
vertical buoyancy frequency, In a thin disk. 


Nt 




(3) 


Fig. [T] maps the vertical shear rate and buoyancy fre¬ 
quency (as well as the gas density) in a fiducial PPD 
model, see 117.1.21 for details. 

Vertical shear is generally weak compared to buoyancy, 
suggesting stability. In fact, without any cooling, the 
Solberg-Hoiland criteria confirms that vertical buoyancy 
i s str ong enough to ensure (axisymmetric) stability, see 
113.41 With radiative cooling, thermal fluctuations decay 
which, combined with pressure equilibration, reduces the 
effective buoyancy. 

How short must tc be for vertical shear to prevail? We 
start with the Richardson number Ri = /{rdz^Y, 

a ratio of buoyant to shear energies. Though not pre¬ 
cisely applicable to our pr oblem, non-rotating sh ear flows 
are stable if Ri > 1/4 (iChandrasekhail HoML also see 
lYoudin fc ShH 120021 : iLee et al.l 120101 f ^ app l icatio ns to 
thin dust l ayers in PPDs). Following lUroinl (j2003fl and 
iTownsencl (jl958l l. we reduce the buoyant energy by the 
ratio of cooling to forcing timescales, tc\rdz^\ < 1 (with 
no reduction in buoyant energy expected when this in¬ 
equality fails). With this reduction, th e Richardso n-like 
criterion for shear instability becomes (jUrpinl [20031 1 


t < (4) 

tc 

If we interpret Eq. [J] as a local criterion, then we see 
that the maximum cooling time which permits instability 
decreases with height as l/|z|; the stabilizing effect of 
vertical buoyancy increases more rapidly away from the 
midplane than the destabilizing effect of vertical shear. 
This finding is relevant to the radiative damping of so 
called ‘surface modes’, which we consider in 1J51 
We are more interested, however, in the global insta¬ 
bility criterion for disturbances at all heights. Our key 
result is the cooling requirement 

tc (5) 

7-1 

for VSI growth that we derive in 115.31 We can simply 
(but non-rigorously) obtain Eq. |S]by evaluating Eq. [T] 
at |z| = 'yH/2. 

We now consider whether, and where, thermal 
timescales in realistic PPDs are short enough to sup¬ 
port VSI growth. Fig. plots the th ermal timescales 
in our fiducial PPD model, see 97.1.21 for details. The 
regions between the red lines satisfy the cooling require¬ 
ment of Eq. 151 Even without a detailed analysis we can 
expect VSI growth between ^ 10 - 150AU. The inner 
disk is complicated by the fact that the optically thick 
midplane has longer cooling times. Thus the cooling cri¬ 
terion fails in the midplane but is satisfied above. This 
case requires the more detailed analysis of m 
Fig. [2] also compares two different radial wavenum¬ 
bers, kx = 10/H and 50/H. In the inner disk, the higher 
wavenumber mode cools faster, due to radiative diffusion. 
Thus in the inner disk, shorter wavelength VSI modes are 
more likely to grow. In the outer disk, the cooling time 
is the same for different kx and z, because optically thin 
cooling is independent of both lengthscale and density. 
Thus the outer barrier to growth, beyond ^ 150 AU in 
this model, is independent of wavelength. This exam¬ 
ple highlights the fact that optically thin cooling sets a 
lower limit to the cooling time; a fact that is sometimes 
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Figure 1. The radial and vertical structure of a ‘minimum mass’ PPD model showing gas density (left, relative to the midplane density 
at 1 AU), vertical shear rate (middle) and buoyancy frequency (right). 
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Figure 2. Estimate of the thermal timescale tc in a fiducial PPD 
model, the Minimum Mass Solar Nebula. Here, /3crit = ^kl/(7~l) 
is a maximum dimensionless cooling timescale below which the VSI 
is expected to operate effectively. This corresponds to the region 
bounded by the red lines. 

obscured when the radiative diffusion approximation is 
made at the outset. 


Cooling times in for k„H=10 


Peril 


Cooling times in for k,<H=50 
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3. GOVERNING EQUATIONS AND DISK MODELS 


An inviscid, non-self-gravitating disk orbiting a cen¬ 
tral star of mass M* obeys the three-dimensional fluid 
equations: 


-f V • (pv) = 0 , 


dp 

m 

dv 1 

--— \-v-'Vv = —VP — 
dt p 

dP 

——h V ■ VP = —yPV ■ V — A, 
dt 


( 6 ) 

(7) 

( 8 ) 


where p is the mass density, v is the velocity field (the 
rotation frequency being 17 = v^/r), P is the pressure, 7 
is the constant adiabatic index, and = —GM,,{r^ + 
13 the gravitational potential of the central star 


with G as the gravitational constant. The cylindrical co¬ 
ordinates (r, (j>, z) are centered on the star. In the energy 
equation (Eq. [S]) the sink term A includes non-adiabatic 
cooling and (negative) heating. The gas temperature T 
obeys the ideal gas equation of state, 

P = -pT, (9) 

where TZ is the gas constant and p is the mean molecular 
weight. 


3.1. Thermal relaxation 

We model radiative effects as thermal relaxation with 
a timescale tc, i.e., 


pTZT-Tcq 

P tc 


( 10 ) 


where Peq = Teq{r, z) is the equilibrium temperature. We 
define the dimensionless cooling time 


/3 = ^Ktc ( 11 ) 

relative to the Keplerian frequency, 17 k = \/ GM^.jr'^. 
We use the terms ‘thermal relaxation’ and ‘cooling’ in¬ 
terchangeably throughout this paper. 

For most of this paper we will take ,9 as a constant 
input parameter, so that we can vary the thermodynamic 
response of the disk from isothermal {jd <C 1 ) to adiabatic 
(/3 1) in a controlled manner. However, in 1J3we will 

consider a more realistic model for jd in PPDs with a dust 
opacity (as in Fig. n. 

3.2. Baroclinic disk equilibria 

The equilibrium disk is steady and axisymmetric with 
density, pressure and rotation profiles p(r, z), P(r, z) and 
17(r, z), respectively. These functions satisfy 


1 dP 


( 12 a) 

pdz 

dz ’ 

1 dP 


( 12 b) 

p dr 

dr 
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A unique solution requires additional assumptions about 
disk structure. We choose 


p(r, 0 ) = pQ{r) = Poo 



T(r, 0 ) = ro(r)=Too , 

P(r, z) = K{r)p^{r, z). 


(13a) 

(13b) 

(13c) 


where tq is a fiducial radius and p and q are the standard 
power-law indices for midplane density and temperature, 
respectively. Without loss of generality we take 9 < 0, as 
is typical in PPDs. The polytropic index P parametrizes 
the vertical stratification with P = 1 (P = 7 ) describing 
vertically isothermal (adiabatic) disks, respectively. The 
ideal gas law requires K = Pq~^TZTq/ p,. 

We further define a modihed sound speed 


c« = (14) 

which in general differs from the isothermal, Ciso = 
Cs/y/T, and adiabatic Cad = Cs y/j/T, sound speeds. 
We also introduce the characteristic scale-height, H = 
Cs(r, 0 )/flK and disk aspect-ratio 

h{r) = — oc (15) 

r 

We are interested in thin disks with h <C 1. 


3.2.1. Density structure 

The equilibrium density field follows from vertical hy¬ 
drostatic equilibrium lEo. I12all . The solution for P 7 ^ 1 
is 


P_ 

Po 


r-i 


1-h 


(r-1) 

*2 


( , ‘ 

y yjl + z^jr"^ 



(16) 


Note that for P > 1 -|- ^,2 there is a disk surface Hg where 
p(r, Hs) = 0; and Hg — y^2/(r — 1)H for h <C 1. 

The vertically isothermal case, P = 1, can be calcu¬ 
lated either as a special case or by taking the limit of 
Eq. [16] as P —1, 


piT,z) 

Poir) 


= exp 


~ exp 


1 


(yrr^ 




for |z| <C r. 


(17a) 

(17b) 


Eq. I17bl is the approximate form of the density field in 
the thin-disk limit. We will primarily focus on vertically 
isothermal disks, the relevant case fo r passively irradi¬ 
ated PPDs (|Chiang fc GoldreichlllQQTH . 


3.2.2. Rotation and vertical shear profiles 

The equilibrium rotation field, D(r, z) f ollow s from the 
density field and centrifugal balance lEa. I12bl) . giving 


D2(r, z) 

D|(r) 


= 1 + 


p + q 
r 


h\r) + 


s 

f 



1 

V^l-hz2/r2 



for all P, where s = q + p(l — P). We also refer to the 
epicyclic frequency n defined via 


2 = ^ ^ ( 4o2t _ 1 ^4^ 

^ ~ r3 ) J.3 Qj. ’ 

where j is the specific angular momentum. 
The vertical shear rate follows from Eq. m 


(3D2 gz 

^ dz hr (1 _p 22 /r 2 )^/^' 


(19) 


( 20 ) 


Eor P = 1, notice that s/P = g is the only disk parameter 
that affects vertical shear. Vertical shear increases lin¬ 
early with height near the midplane, and | 9 zD 2 | is max¬ 
imized at min(r/-\/ 2 , Hg), which is expected to limit the 
maximum VSI growth rate. 

A more general expression for vertical shear holds for 
vertical polytropes (which satisfy Eq. I13cl) , with no as¬ 
sumed radial structure: 

din p dK g^dlnK 

’ Ih =^ —17 =-T—’ 

where pz = dzPjp. 


3.3. Entropy gradients and vertical buoyancy 

The gradients of specific entropy, S = Cp InP^P/p, in 
our disk models are 


or ijr 


7 


- 1 


ainp 

dz 


dlnp 


dr 

Cppz r - 7 

ci 7 


( 22 a) 

( 22 b) 


where Cp is the heat capacity at constant pressure. 
The vertical buoyancy frequency Nz is 


Nl = -Cfi^gz^ = 


dz 


cl 7 


(23) 


We only consider convectively stable disks with with P < 
7 and Nl > 0 . 


3.4. Stability without cooling 

In the absence of cooling, axisymmetric stability is en¬ 
sured if both Solberg-Hoiland criteria are satisfied: 


dP 

dz 


pCp 

_ 

dr dz 


VP 

• VS' > 0, 

(24a) 

dz 

dr)^ 

(24b) 


(|Tassoullll978f l. Eor rotationally-supported, thin disks 
the first criterion is easily satisfied. Thus, we consider 
the seco nd criterion. Without loss of generality (since 
Eq. I24bl is even in z), we consider z > 0 so that dzP < 0. 
With drj"^ — and using Eq. I^Ijfor dzj^ we find that 


7 — P > 


P 



s( 7 -P) 


dlnp 

dlnr 


(25) 


implies stability. For typical model parameters, the 
right hand side is 0{h^) < 10“^. The left hand side 
is positive and order unity in our disk models (with 
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7 — r > 0.1) implying strong stability to convection, see 
Eq. [23l Th is stable stratification is e xpected for irradi¬ 
ated disks (|Chiang fc GoldreichlllQQTH . Thus, adiabatic 
disturbances in a standard, irradiated disk are stable to 
a disk’s vertical shear, explaining why the VSI requires 
cooling. 


4. LINEAR PROBLEM 

This section presents the two sets of equations we use 
to study the linear development of the VSI. Both sets are 
radially local and vertically global with a finite cooling 
time. The first set, presented in M.ll is more general and 
used for numerical calculations. The second set, in M.21 
makes additional approximations about disk structure 
and wave frequency. This simplified set is used to obtain 
analytic results. 


4.1. Radially local approximation 

We consider axisymmetric perturbations to the above 
disk equilibria. The growth of the VSI is strongest for 
short radial wav elengths, as compared to the disk ra¬ 
dius (|N13l : IBLISH . We thus perform a standard two step 
process to obtain linear ized equations in the radial ly lo¬ 
cal approximation fe.g.. lGoldreich &: Schuber^ll967L who 
also consider vertically localized perturbations). 

We first expand all fluid variables into the equilibrium 
value plus an Eulerian perturbation denoted by i5, e.g. 
5p for the perturbed density, and drop all non-linear 
perturbations. Second, we perform a Taylor expansion 
about a fiducial radius rp with the local radial coordi¬ 
nate X = r — rg, keeping only the leading order terms in 
x/xq. We also relabel (trivially for axisymmetric pertur¬ 
bations) the azimuthal direction as y. 

Perturbations take the form of a radial plane wave with 
arbitrary vertical dependence, e.g. 

SpSp(z) exp (ikxX — ivt), (26) 


where v = w -|- icr is a (generally) complex frequency, 
with w and a being the real frequency and growth rate, 
respectively, and kx is a real radial wavenumber. We 
take kx > 0 without loss of generality, assume kxXg ^ 1 
and neglect curvature terms. Henceforth all unperturbed 
fluid variables, including gradients such as refer to 
equilibrium values at (ro,z). 

We further define the perturbation variables W = 
6P/p and Q = c^^SpIp. With this procedure the lin¬ 
earized system of equations are 


. Q -7, A 

IV— = ikxOVx + 

d 


(? In p „ i9 In p „ 

-j - 1 - TT^dvz + C q ^ 5vx, (27a) 

dz az or 


\v5vx = ikxW — 2Vl5vy — 

^ ^ p dr cl 

dfl K? 

IVOVy = rg — OVz + TX^OVx, 


dz 


2H 


, dW dlnp 

iv5vz = — -k (W-Q), 

ivW = (^^kx6vx + 


dSvz 

dz 


W- 


Q 


+ C-^dVx- 
p dr 


(27b) 

(27c) 

(27d) 

(27e) 


This is a set of ordinary differential equations (ODEs) in 
z. Solutions to these are presented in )j6] and an The 
coefficient C = 1 is introduced simply to label terms with 
explicit radial gradients of the equilibrium state, which 
again are evaluated at (ro,z). For clarity, hereafter we 
drop the subscript 0. We will consider the effects of ig¬ 
noring these radial gradients below. 

4.2. Reduced equation for vertically isothermal disks 

Our simplified model starts with Eqs. [23 and makes 
the following additional simplifications: 

1. We set r = 1, focusing on vertically isothermal 
disks. 


2. We set C = 0, neglecting terms with an explicit de¬ 
pendence on the radial structure of the equilibrium 
disk. Vertical shear, which implicitly depends on 
the radial temperature gradient, is retained. This 
fully-radially-local approximation is also m ade in 
the ‘vertically global shearing box ’ of IMP 14 

3. We make the low frequency approximation, as¬ 
suming \ v‘^\ <C K^, Si milar to th e incom¬ 

pressible (jGSBTf ) or anelastic (|N13l : IBL15I 1 approx¬ 
imations, the low frequency approximation filters 
acoustic waves in favor of inertial-gravity waves 
(iLubow fc PringldfTM^ . 

4. We make the Keplerian approximation, setting 
H —>■ Hk and k —>■ Hk, but retaining the verti¬ 
cal dependence in the crucial vertical shear term, 

dz^. 


5. We consider thin disks with h <C 1 and use the 
Gaussian approximation, Eq. I17bl for the equilib¬ 
rium density field. 


In terms of the dimensionless variables 

z = z/i7, k = kxH, V = Lu + id = v/Hk, (28) 

where ui and d are real, the above approximations lead 
to a single second order ODE, 


Sv" — zA6v'^ + {B — Cz^)5vz = 0, 

(29) 

where ' denotes d/dz and 


A = 1 + ihqit, 

(30a) 

B = v'^ (^X + - (x + , 

(30b) 

:S 

1 

1 

III 

(30c) 

where 


1 — iv/3 
^ 1 — iv/3"/ 

(31) 

The derivation of Eq. [23 is detailed in Appendix [^ 

In Appendix |B] we discuss the limits of these approx¬ 
imations. In particular we point out that the fully- 
radially-local approximation is only valid for short cool¬ 
ing times, d < 0(1), which is the regime in which the VSI 
operates, as demonstrated below. On the other hand, for 


d > 0(1), this approximation can introduce artificial in¬ 
stability due to the n on-s elf-consistent neglect of global 
radial gradients (see EU. 
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5. ANALYTIC RESULTS WITH FINITE COOLING TIMES 

Previous analytic studies of the VSI have largely fo¬ 
cused on isothermal perturbations, with infinitely rapid 
cooling, as discussed in the introduction. We further ex¬ 
plore this /3 = 0 limit in Appendix [C] both to further 
develop intuition for this idealized case and to establish 
a connection with previous works. However our main 
interest is the effect of finite cooling times, which we ex¬ 
plore analytically in this section. 

In il5.ll we show that even an infinitesimal increase in 
the cooling time, starting from /3 = 0, slows the growth 
rate of the VSI. In il5.2l we derive exact solutions to the 
simplified VSI model developed above (Eq. [^n]) . In il5.3l 
we reach our main result, the maximum cooling time 
above which VSI growth is suppressed. 


5.1. Introducing a small but finite cooling time 

We are interested in finite thermal relaxation 
timescales /3 > 0, but it is instructive to first ask the 
more analytically tractable question: how do the eigen- 
frequencies and eigenfunctions change when we change 
(3 from zero to a small but finite value? For sufficiently 
small P we expect the solution to only differ slightly from 
a case with /3 = 0. We thus perturb a solution for /3 = 0 
to see the effect of finite cooling. 

For definiteness, let us consider the simple solution 


/3 = 0, 5vz = 1, 


1 -|- ihqk 
1 + F ’ 


(32) 


which solves Eq. [29] since 5vz =constant and x = 1 so 
that B = C = 0. This is the lowest order VSI mode, the 
‘fundamental corrugation mode’, where the vertical ve¬ 
locity is constant throughout the disk column. Hereafter 
we shall simply call it the fundamental mode. 

The fundamental mo de has been observed to dominate 
numerical simulations ([Stoll fc Klevll20l4 IN13I1 . and are 
the ones we hnd to typically dominate in numerical cal¬ 
culations with increasing /? f ( 16.211 for moderate values of 
k, as well as in PPDs with a realistic estimate for ther¬ 
mal timescales ((JT]). We will see later that consideration 
of low order modes in fact provides a useful way to char¬ 
acterize the effect of increasing the thermal timescale on 
the VSI. 

We linearize Eq. [53 about the above solution for /3 = 0 
and write 


/3 —>■ 0 -I- <5/3, Svz —t Svz + Svzi, v ^ V + Sv, (33) 
which implies 

X^l + 6x = l + ivi'y-l)6p, (34) 

with Svz and v given by Eq. 1321 We may then seek 

5vzi = az^ + b, (35) 

where a, b are constants. 

We insert Eq. |33| into Eq. |29l keeping only first order 
terms, and solve for Sv using the above expressions for 
6x and v. We find imaginary part of Sv is 

(7 — 1 ) (k"^ — 2h'^q^k'^ — lPq^\ 

Sa = ----- 2 ^ Sp. (36) 

2 (1 + K^q^k^] (1 -I- ) 


Since /i <C 1, introducing finite cooling SP > 0 implies 
Sa < 0, i.e. stabilization, since 7 > 1. 

A finite cooling time allow buoyancy forces to stabilize 
vertical motions in sub-adiabatically stratified disks. The 
dependence in Svzi makes sense because it becomes 
significant at large |z|, i.e. away from the midplane where 
the effect of buoyancy first appears as p is increased from 
zero (since Nf oc z"^ for a thin disk, see Eq. (231) . 

5.2. Explicit solutions and dispersion relation 
We now solve Eq. [53 explicitly. We first write 

duz(£) =g(z)exp , (37) 

where A is a constant to be chosen for convenience. In¬ 
serting Eq. |37| into Eq. [551 gives 

0 = g'' -z {A- 2A) g' + {B + X)g+ (A^ -XA-C) z^g. 

(38) 

We choose A to make the coefficient of z'^g vanish, and 
impose the vertical kinetic energy density p\Svz\'^ oc 
Igp exp (Re A — 1/2) to remain finite as |z| —>■ 00 . 
Then assuming g[z) is a polynomial, we require 


Re A < i, (39) 

which amounts to choosing 

A = i (a - \JA^ -L 4c'j . (40) 

Eq. [29l becomes 

0 = g” — z {A — 2X) g' + {B + A) g. (41) 
We seek polynomial solutions 

M 

g{z) = ^ b^z^, (42) 

m—0 

which requires 

B{v) + X{v) =M[A- 2A(b)]. (43) 

For ease of analysis, we re-arrange Eq. |43| and square it 
to obtain a polynomial in v, 

6 

0 = ^c/t)', (44) 

1=0 


where the coefficients ci are given in Annendix ID.ll The 
dispersion relation v = v{k) is complicated and gener¬ 
ally requires a numerical solution. However, simple re¬ 
sults may be obtained in limiting cases which we consider 
below. 


5.3. The critical cooling time derived 

Here we estimate the maximum thermal relaxation 
time scale Pc that allows growth of the VSI. In Appendix 
ID. 21 we derive the following relation between Pc and the 
wave frequency Wc, assuming marginal stability to the 
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VSI: 

(^LUckj = (^ujck^ + M{M + 1) ^1 — , (45a) 

=^i"f - 1)(1 + 2M)^ (^tockj (45b) 

- 2M(M + 1). 

These equations consider k'^ 1 and assumes \dik\ is 

0(1) and Pc "C 1. Note that Eq. I45al reduces to the 
dispersion relation for low-frequency inertial waves in the 
absence of vertical shear (see Appendix I C. 2.1 1 for details). 

We are most interested in the longest cooli ng t ime 
which allows growth for any M. In Appendix ID.3I we 
find that if \hqk\ < 1, so that marginal stability ex¬ 
ists, then Pc decreases with increasing M. The VSI 
criterion is then set by Pc for the fundamental mode 
(M = 0), for which the exact solution to Eq. I45al — 

I45bl is uJck = Pdj — i)/hq = —1. We thus find the cool¬ 
ing criterion for VSI growth is 

P < Pcrit = (46) 

7-1 

The thin disk approximation, h <C I, indi cates that 
/3crit <C 1, as assumed to obtain Eq. I45al — libbl We 
heuristically explain Pait in ^ and numerically test its 
validity in ^6.31 

6. NUMERICAL RESULTS FOR FIXED COOLING TIMES 
This section investigates the linear growth of the VSI 
with fixed cooling times, that are independent of 

height. We solve the radially local model of Eos. E71 with 
C = 1, i.e. retaining radial gradients of the background 
disk structure. The equilibrium disk structure given is in 
>13.21 This numerical study does not make the simplifying 
approximations used in ^ 

We solve the linear problem by expanding the per¬ 
turbations in Chebyshev polynomials T; up to 1 = 512 
and discretizing the equations on a grid with z C 
[—•^max) Zjnax]- Our Standard boundary conditions at the 
vertical boundaries is a free surface, 

AP = SP + ^ - VP = 0 at z = izmax, (47) 

where ^ is the Lagrangian displacement with meridional 
components ^x,z = ^Svx,z/'v. In some cases we impose a 
rigid boundary with (5uz(±Zniax) = 0. 

The above discretization procedure converts the linear 
system of differential equations to a set of algebraic equa¬ 
tions, which we solve using LAPACK matrix routine^ 
Unless otherwise stated, we adopt fiducial disk parame¬ 
ters ( 7 , r) = (1.4,1.011), (p, q, h) = (—1.5, —1,0.05), and 
vertical domain size Zmax = 5i7. This setup is effectively 
vertically isothermal since c^Zmax) — 0.9Cg(0)EI In Fig. 
[3] we plot the corresponding equilibrium rotation profile, 
which shows that Vl decreases slightly away from the mid¬ 
plane. From Eq. 021 we expect that P < /3crit = 0.125 is 
needed for VSI growth in our fiducial model. 

^ Available at http://www.netlib.org/lapack/ 

® While this mod el has a zero density surface at z = ±Hs ~ 
±16.4// (see Eg. II 8 II . the surface is outside the numerical domain. 



z/H 


Figure 3. Equilibrium rotation profile fl(z), normalized by its 
mid-plane value, for the fiducial disk model with E = 1.011 and 
(p,g, h) = (-1.5,-1,0.05). 



-a)/(hfl|^) 

Figure 4. Growth rate, a, vs. oscillation frequency, to, for VSI 
modes with wavenumber kx = k/H = 10/H in the fiducial disk 
model with rapid cooling, /3 = 10 ~®. We investigate the effect of 
vertical box size, from Zmax = OH [blue triangles) to 7H {green 
x’s) with a free surface boundary (AP = 0 ). For Zmax = 5// 
(black dots) we also compare free surface {black dots) and rigid 
{red plusses), finding little difference. The fundamental mode is 
marked by ‘f’. Lines show the analytic predictions for isothermal 
perturbations (/3 = 0 , dashed) and for fi = 10 “® {dotted). See text 
for discussion. 

6.1. Rapid thermal relaxation 

We first calculate the VSI in a disk with rapid thermal 
relaxation by setting P = 10“^. Since P <C /3crit, this case 
gives similar result s to the well studied case of isothermal 
perturbations (e.g. IN13I : IMP14 lBL15f l 0 This case allows 
us to explore and test our finite cooling time model in a 
familiar context. 

6.1.1. Case study: fc = 10 

Fig. 0]maps the growth rates a and wave frequencies w 
of unstable modes for k = 10. Each symbol corresponds 
to a different order mode, i.e. a different number of ver¬ 
tical node crossings. Lower order modes have smaller 
|a;| for fixed radial wavenu mber, as expected for inertial- 
gravity waves (see 8C.2.1I1 . 

We compare our numerical results to the ana lytic dis¬ 
persion relations for /3 = 0 and 10“^, from Eos. IC16l and 
[44l respectively. The analytic models also have discrete 
modes, which lie on the continuous curves that are plot¬ 
ted for clarity in Fig. 01 While the analytic models have 
many simplifications, they crucially lack an artificial ver¬ 
tical surface (as they assume an infinite vertical domain). 

^ We confirmed that smaller /3 values gave similar results. 
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co= -2.04hnK, o=0.45hnK, (3=0.001, k,H=10.0 



z/H 

Figure 5. Pseudo-enthalpy perturbation IV (top) and vertical 
velocity perturbation Sv^ of the fundamental VSI with fc = fO in 
the fiducial disk with rapid thermal relaxation /3 = 10~®. The real 
(imaginary) parts of W and Sv^ are plotted as solid (dotted) lines. 




0 2 4 6 8 10 12 14 
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Figure 7. Same as Fig. [4] but for k = 30. Examples of surface 
modes are marked by ‘s’. 

rigid boundaries (red plusses). 

Rather, the decline in growth rate for high order mod¬ 
els is due to the vertical box size, as seen by compar¬ 
ing the z-catajH = 3, 5 and 7 cases in Fig. (H Increas¬ 
ing Zniax give better agreement with the analytic theory, 
which have Zmax —t oo. Larger boxes include regions of 
larger vertical shear (see Eqs. El and EOl) . which higher 
order modes can tap to give larger growth rates. For 
the 2max = 7iL case, we begin to see another branch of 
modes with the highest growth rates at a; ^ —G/iIIk- 
This branch contains the ‘surface modes’ described be¬ 
low. 

The need for vertically extended domains to capture 
the largest VSI growth rates is problematic, especially 
for hydrodynamic simulations. This complication is mit¬ 
igated by at least two factors. First, we show below that 
with slower cooling, the growth of higher order modes 
is preferentially damped, consistent with the analysis in 
>15.31 Second, the fastest growing modes may not dom¬ 
inate transport when they operate in very low density 
surface layers, i.e. at many H. Quantifying which modes 
will contribute most to non-linear transport is an impor¬ 
tant issue, but beyond the scope of this work. 

Fig. 0 shows eigenfunctions W = 5P/p and 5vz for 
the lowest order fundamental mode. (In this and similar 
plots below, we normalize eigenfunctions such that its 
maximum amplitude is unity with vanishing imaginary 
part at the lower boundary.) Fig. [5] maps the pressure 
perturbation and meridional flow, scaled by y/p to reflect 
the contribution to kinetic energy. Notice the stretched x 
axis. Radial velocities are in fact typically much smaller 
than vertical velocities, as e xpected for a vertically elon¬ 
gated, anelastic mode (IN13I1 . Most of the kinetic energy 
is contained within ~ 2H of the midplane due to the 
density stratification. 


Figure 6. Visualization of the fundamental VSI mode in Fig. [S] 
For this plot the mode has been transformed back to real space. 
Thus, the color scale corresponds to the real pressure perturbation 
scaled by its maximum value. Arrows show the real meridional 
flow multiplied by y/p. The horizontal axis has been stretched for 
clarity. 

For our fiducial case (black dots), the lower order and 
lower frequency modes closely follow the analytic pre¬ 
diction with growth rates increasing with |w|. However, 
for |a;| > 4 /iHk the growth rate declines with |a;|. This 
break from the analytic prediction is not due the choice 
of boundary condition, as we demonstrate by considering 


6.1.2. Surface modes of the VSI 
The surface modes mentioned above are more promi¬ 
nent for larger wavenumbers. Fig. [7] maps the fc = 30 
eigenvalues, with the surface modes labeled. Surface 
modes strongly depend on the location of the imposed 
vertical boundary, and disagree with the analytic mod¬ 
els, which lack an imposed surface. We thus discount the 
physical significance of surface modes for our model, as 
we explain further below. 

Surface modes are a well known fe at ure of the VSI in 
finite vertical domains (|NI3t IMPI4D . IBL15I show that 
surface modes arise whether the surface is imposed (as 
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co= -1 .SSh^K- 0=1 ■ 48 hnK, ( 3 = 0 . 001 , k,H= 30 



z/H 

Figure 8. Example of a ‘surface mode’ in the fiducial disk model. 

in our models) or natural (for disk models with a zero- 
density surface). BL15 mention the interface between a 
disk’s interior and corona as a physical boundary that 
can trigger surface modes. Since vertically isothermal 
disk models lack a surface, modes which depend on the 
adopted value of Zmax are not physically meaningful. 

The possibility of artificial surface modes seeding 
growth in numerica l sim ulations of the VSI merits fur¬ 
ther study. Indeed IN13I find that the initial growth of 
perturbations primarily occurs near the vertical bound¬ 
aries. The motion in surface modes is indeed concen¬ 
trated near the surface, as shown in Fig. [SJ where the 
density is low. Thus, their contribution to transport 
might (by themselves) be weak, following the arguments 
in ilb.l.ll Moreover surface modes, like all modes with 
large wavenumbers, are more prone to viscous damping, 
as discussed in ^8.11 

While the relatively large growth rates of surface 
modes is tantalizing, we dismiss them in disk models that 
lack a physical surface. 

6.2. Slower thermal relaxation 

We now consider the effect of longer cooling timescales 
by gradually increasing /3. We expect VSI growth for 
13 < /3crit = 0.125. We also expect that higher order 
modes will damp at even lower (3 values, as argued in 

Fig. |9] largely confirms these expectations by plotting 
eigenvalues for fc = 10 and (3 from 0.01 to 0.1. The 
analytic results from Eq. |44] are now plotted as different 
curves for each mode order M, with f3 varying along each 
curve. We see the standard increase in frequency, |w|, 
with mode order. 

As expected, the higher order modes are preferentially 
damped. For (3 > 0.03 the fundamental [M = 0) mode is 
the fastest growing. For (3 = 0.1, only the M = 0 mode 
grows. This growth is slow since /3 is near Pent- 

Fig. |9] also shows how cooling times affect the de¬ 
pendence on Zmax, the size of the vertical domain. For 
all cooling times, the fundamental mode depends only 
weakly on Zmax, and there is good agreement with an¬ 
alytic values. This convergence is reassuring given the 
importance of the fundamental mode at longer cooling 
times. For higher order modes and short cooling times, 
however, the eigenvalues vary strongly with Zmax- The 
disagreement with analytic theory is strongest for the 
smallest domain, with Zniax = 3iJ. For /3 > 0.03, the 
ZmaxIH = 5, 7 and the analytic results are well converged 
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Figure 9. Unstable modes with fc = 10 and thermal relaxation 
timescales /3 = 0.01 (black dots), /? = 0.03 (red crosses), /3 = 0.05 
(blue triangles) and /3 = 0.1 (green crosses); for different vertical 
domain sizes Zmax = 3H (top), Zmax = SV (middle, the fiducial 
set up) and ^max = 7H (bottom). Lines are computed from Eq. 
1441 for M = 0 (solid, fundamental mode) and M = 1, 2, 3 (dotted, 
dashed, dash-dot, respectively). Along each line, j3 increases con¬ 
tinuously from 0.01 to 0.1 from top to bottom, and squares mark 
corresponding 0 values with numerical results. 

(aided by the fact that only M < 5 modes grow). 

Figures [TU] and [TT] illustrate the eigenfunction of the 
fundamental mode with fc = 10 and /3 = 0.1. Com¬ 
pared with the small /3 case in Figs. [S] and [HI the eigen- 
fuctions show a more complex variation of phase with 
height. This dependence yields the ‘tilted’ appearance of 
pressure field in Fig. [TT] Physically, the increased role 
of buoyancy, which increases with height, explains why 
larger /3 values produce more complex vertical structure. 

Fig. IT^ shows how cooling affects higher wavenumber 
modes, specifically k = 30. We only show the larger 
domains with Zmax = 5iT, 7H. These cases again show 
good convergence for /3 > 0.03, which the smaller domain 
with Zmax = 3iT lacks. 

With slower cooling, the wave frequency shows a more 
complex dependence on mode order, both analytically 
and numerically. For /3 = 0.1, the fundamental mode 
(which lies on the M = 0 curve) no longer has the small¬ 
est I a; I value. 

Moreover, the fundamental mode is no longer the 



































10 


M-K. Lin, A. N. Youdin 



-4 -2 0 2 4 

z/H 


Figure 10. Same as Fig. [5] but with a thermal relaxation 
timescale /3 = 0.1. 



Figure 11. Same as Fig. [6] but with a thermal relaxation 
timescale /3 = 0.1. 


fastest growing mode for /3 > 0.03, or even for /3 = 0.1. 
This complication is not actually surprising since \ hqk\ = 
1.5 is no longer less than unity, as required in the analytic 
derivation of 115.31 Though inconvenient, given that the 
fundamental mode is easy to identify and the most nu¬ 
merically converged, this complication is not ultimately 
significant for the operability of the VSl. 

Fig. [T^lalso demonstrates that artificial surface modes 
are damped for P > 0.03. We are thus confident that ar¬ 
tificial surface modes should not affect our determination 
of the critical cooling time. 
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Figure 12. Same as Fig. but for fc = 30 and domain sizes 
Zmmx = 5H {top) and 7H {bottom). Examples of surface modes 
are marked with ‘s’. 



Dimensionless cooling time p 


Figure 13. Maximum VSI growth rates in the fiducial disk model 
as a function of the thermal relaxation timescale 0, for Zmax = 5H 
(top) and Zmax = 7H (bottom). The ve rtica l line is the critical 
thermal timescale /9crit obtained from Eq. 1461 

6.3. Critical thermal relaxation timescale 

Having explored the behavior of VSI modes with cool¬ 
ing, we turn to the numerical validity of the analytic cool¬ 
ing criterion for vertically isothermal disks, /3 < /3crit = 

Fig. [T51 shows how VSI growth rates vary with cooling 
time in our fiducial model with /3crit = 0.125. Curves are 
for a hxed horizontal wavenumber, and show the maxi- 
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mum growth rate for all vertical mode orders. The dis¬ 
continuity in some curves occurs when the fastest growth 
switches to a different mode order. 

For k = 5,10, the growth rate drops to zero at the 
expected /3 = /3crit- For longer wavelength modes, with 
fc = 1 , growth persists to slightly longer cooling times. 
This difference is not surprising since our analytic deriva¬ 
tion assumed ^ 1 , but the change is quantitatively 
minor. 

For shorter wavelengths, with k > 30, growth persists 
for P significantly larger than ,5crit- This tail of growth is 
partly explained by the breaking of the k < \/\qh\ = 20 
approximation, used in the analytic derivation. Despite 
the lack of a clear stability boundary at high fc, the 
Pcvit threshold remains useful, since growth rates drop 
to < 10 % of their maximum value at /3crit and continue 
to fall for larger /3. Moreover, we expect that longer 
wavelength modes with fc < 20 are more significant for 
disk transport, see ES 

Comparing the Zmax/^f = 5,7 cases in Fig. [131 we 
see that the location of the vertical boundary has little 
effect on the critical cooling time. This agreement occurs 
despite the fact that peak growth rates differ as ,5 —>■ 0. 
We are thus confident that boundary effects, including 
surface modes, do not affect our analysis of the critical 
cooling time. 

Our res ults a gree with the vertically isothermal simu¬ 
lations of IN13I which used the same disk parameters as 
our fiducial model. The expected ,5crit = 0.125 is con- 
sistent with the simulations shown in their Fig. 12 . IN13I 
found nonlinear VSI growth for /3 = 0.06 < Petit but no 
growth for /3 = 0.6 > /3crit0 

Fig. 1141 confirms that the numerically determined Petit 
shows the expected scaling with disk parameters h, q, and 
7 . For this test we fix fc = 10 (an appropriate value for 
all the reasons discussed above) and measure the small¬ 
est P value where growth vanishes. The agreement with 
the analytic scaling of Eq. 0S] is quite good, confirming 
the applicability of our critical cooling time in vertically 
isothermal disks. 

6.4. Vertically non-isothermal disk model 

We now generalize our critical cooling time by con¬ 
sidering disks which are not vertically isothermal. We 
rely on physical arguments for this generalization. In 
disks with weaker vertical stratification, i.e. F > 1 for 
fixed 7 , we expect Petit to be larger. We thus transform 

1 /( 7 - 1 )^ 1 /( 7 -F). 

We also expect petit to scale with the vertical shear. 
In general the vertical shear rate dzVl oc s, the radial 
entropy gradient, see Eq. |20l We thus transform q —^ 

s = 9 + p(i-r). 

Our estimate for the generalized cooling time crite¬ 
rion — valid for both vertically isothermal and non- 
isothermal disks (cf. Eq. H5)) — is thus 

Petit ^ /5crit,gen — p; ■ 

7-1 

We test /3crit,gen using a disk model with {p,q,h) = 
(—0.5,—1,0.05) and ( 7 ,r) = (1.4,1.3). For this setup, 

^ Since the dimensionless cooling time T^eiax in lN13l is normalized 
to the orbital period, we convert P = f^K^orb^relax = 27rTreiax- 



|q| 
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h 

Figure 14. Dependence of the upper limit to the thermal relax¬ 
ation timescale /3crit foi" fhe fundamental VSI mode on disk pa¬ 
rameters. The fiducial setup is ( 7 , F = (1.4,1.011) and (p,q,h) = 
(—1.5,—1, 0.05). Top: varying g G [—1.2,—0.2]; middle: varying 
7 G [1.2, 2.0]; bottom: varying h G [0.02,0.1]. The perturbation 
wavenumber is k = 10 . 

s = —0.85 and /3crit,gen = 0.425, which is ^ 3 times 
larger than our fiducial, vertically isothermal disk with 
^crit = 0.125. 

We set the vertical domain size just inside the physical, 
zero-density disk surface, 2 „iax = 0.9977s ~ 2.677, see Eq. 
m At this surface, the vertical buoyancy frequency, TV^, 
diverges as \z\ —^ Hg. While we might expect patholog¬ 
ical behavior for this model, we fortunately do not find 
it. 

In Eig. [T5]we plot the mode diagram for k = 30 and 
several values of p. This plot is qualitatively similar to 
the vertically isothermal case, Eig. 0 As before, larger 
P values rapidly stabilize higher order modes. Eor suf¬ 
ficiently large P only the fundamental mode is unsta¬ 
ble, with a growth rate that decreases as P approaches 

^crit.ge n- 

Fig. [in] shows how maximum VSI growth rates vary 
with P in the vertically non-isothermal model. The qual¬ 
itative behavior is similar to Fig. [T3| for the vertically 
isothermal disk. Eor k = 0(10), growth rates drop to 
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7.1.1. Radiative diffusion and Newtonian cooling 

The thermal relaxation of a perturbation depends on 
the relative sizes of the perturbation’s lengthscale, I, and 
the photon mean-free-path, 

= — (49) 

KdP 

where Kd is the opacity, specifically the dust opacity ap¬ 
propriate for cold PPDs. 

In the optically thick limit, 1 Irad, radiative diffusion 
smooths out thermal perturbations. In this regime, the 
linearized cooling function is 


Figure 15. Unstable modes in a vertically non-isothermal disk 
with /c = 30 and a range of thermal relaxation timescales. 



Dimensionless cooling time p 

Figure 16. Maximum VSI growth rates in the vertically non- 
isothermal disk model as a function of the thermal relaxation 
timescale tc = The ve rtica l dashed line is the critical ther¬ 

mal timescale given by Eq. \M A conjectured critical thermal 
timescale accounting for a n on-isothermal background, given by 
/3crit,gen = 0.425 (see Eq. 1 4811 . is also shown as the vertical dotted 
line. 

zero at the predicted ^Scrit.gen- As in the vertically isother¬ 
mal case, modes with both lower and higher wavenum¬ 
bers exhibit slow growth rates beyond the critical cooling 
time. As argued in H6.31 the critical cooling time is useful 
despite not being an absolute stability limit. 

The generalized cooling time criterion in Eq. 05] thus 
appears to be valid. If a disk’s vertical structure is not 
characterized by a single polytropic index, P, we expect 
that a density weighted average of /3crit,gen should be a 
good approximation. This expectation remains to be 
tested. 

7. THE VSI WITH REALISTIC PROTOPLANETARY DISK 
COOLING 

To this point we have treated the cooling time, /3, as a 
free parameter that, moreover, is independent of height. 
To understand the applicability of the VSI to PPDs, we 
now consider /3 values that are consistent with PPD mod¬ 
els. Section [70] develops our PPD cooling law, which de¬ 
pends on disk radius, disk height and perturbation wave¬ 
length. We compare these PPD (3 values to /3crit in ^7.21 
While this comparison is instructive, it is not complete 
since /3crit was derived assuming a vertically constant /3. 
Thus in 117.31 we analyze the linear growth of VSI with 
our PPD f3 values. 


7.1. Cooling model 


Mdiff = ■ (A:,adV5T), (50) 

pTCy 

where fcrad is the radiative conduction coefficient defined 
below, and Cy is the specific heat capacity at constant 
volume. 

Since the VSI is characterized by vertically-elongated, 
radially-narrow disturbances, we retain only the radial 
derivatives of the perturbations in Eq. 1501 
radially local approximation, we have 

St 

dAdiff ~ -pklP—, 

where the radiative diffusion coefficient is 

^ fcrad _ IQctsT^ 

^ pCy iKdP^Cy ’ 

and fJs is the Stefan-Boltzmann constant, 
relaxation timescale (defined in Eq. flUl) 
diffusion is thus 

rjkl' 

In the optically thin regime, I <C Zrad 
taxation operates by ‘Newtonian cooling’, 
time is independent of I and inversely proportional to 
the opacity, tthin oc f/nd- Specifically 

ithin = (54) 

67] 

and Ahin does not depend on p (because our adopted Kd 
depends on T only, see below). 

Our general cooling time for all perturbations, 

ic = ^thin 4“ Idiff j (55) 

is a simple prescription to smoothly match the optically 
thick and thin limits. 


Thus, in the 

(51) 

(52) 

The thermal 
for radiative 

(53) 

, thermal re- 
The cooling 


7.1.2. PPD cooling times 

In a vertically isothermal disk with surface density S = 
PoV^H, the dimensionless thermal relaxation time /3 
becomes 


I3{r,z;k) = 


r 1 pHzY 

?7p2 |_3 k^E 2 27rfc2 ’ 


(56) 


where p = p/po- The first and second terms in square 
brackets represent the optically thin and thick cooling 
regimes, respectively. 
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z/H 

Figure 17. Thermal relaxation timescales in the fiducial MMSN 
at r = 50AU and r = 5AU for k = 30 (solid lines). The corre¬ 
sponding horizontal dashed lines are the critical thermal relaxation 
timescales derived in linear theory. 


We adopt the Minimum Mass Solar Nebula (MMSN) 
disk model of iChiang fc Youdinl (120101) which specifies 


E = 2200r“u/^gcm-^ (57a) 

T = 120rXu'^^K, (57b) 

with TAu = r/AU. This model has q = —3/7, p = 

—39/14 and h = 0.022r^(J, using p = 2.33 and the above 
relation between E and po rx r^. 

For the dust opacity we adopt the iBell fc Lhil (|1994[) 
law with Kd oc. T^, giving 


Kd = 2.88kdr’Au^’^cm^ g \ (58) 


where the opacity normalization factor, scales with 
the ratio of small dust to gas; Kd = 1 is our fiducial case. 

For this disk and opacity model, the cooling time be¬ 
comes 


/3(r, z; k) =3.9 x 10 


1 -h 1.9 X Wklp'^{z)r^^'''‘ k 


-33/71-2 


(59) 


Perturbations are in the optically thin regime for 

k > 2.5 X 10^Kdp(z)r^u^^^^. (60) 

In the inner disk, e.g. tau 1; only extremely small- 
scale perturbations are in the optically-thin regime near 
the midplane. However, in the outer disk, e.g. tau 
10, more moderate wavenumbers fc ^ 10 experience 
optically-thin cooling. 


7.2. /3 vs. Pent in PPDs 

The simplest way to estimate whether the VSI can op¬ 
erate at a given radius in a PPD is to compare the cooling 
time (Eq. [59]) to its critical value (Eq. |4l|) , evaluated for 
the MMSN as 

Pent = 0.0244/J. (61) 

This comparison is strictly valid only for optically thin 
cooling which is independent of height, as assumed in the 
derivation of Pent- This condition typically holds in the 
outer disk. Fig. |T7| shows that a fc = 30 perturbation 
at 50AU has vertically constant p. Furthermore, since 
P < Pent, the disturbance would be unstable. 


For optically thick cooling, P varies with height, com¬ 
plicating the comparison with Pent- At 5AU, Fig. |T7| 
shows that cooling times are long in the midplane with 
P > Pent- However, since P < P ent a way from the mid¬ 
plane, we require the analysis of >17.31 to determine if the 
VSI can grow. (That calculation will show that growth is 
strongly suppressed for this case.) We proceed with the 
awareness that optically thick regions in the inner disk 
are a complication, but that we can reasonably expect 
VSI growth if /3 < Pent at all heights. 

In Fig. IlHl we compare P to Pent across a range of disk 
radii for different heights, opacity values and wavenum¬ 
bers. For optically thin cooling in the outer disk, curves 
for different wavenumbers overlap, as expected from Eq. 
(591 Since the (optically thin) slope of P is steeper than 
for /3crit, VSI growth can be suppressed at large radii (for 
the chosen opacity law). This effect is seen for = 1 
in the central panels of Fig. [181 where VSI is damped 
outside ^ 150 AU. 

We find that the most important factor for VSI growth 
is the opacity. With a smaller opacity, kd = 0.1, growth 
is suppressed at all radii, as shown in the top panels of 
Fig. [iHl Since optically thin cooling is too slow in this 
case, optically thick cooling — above the floor set by 
optically thin cooling — is also too slow. 

Larger opacities make optically thin cooling much 
faster than Pent, as shown in the top panels of Fig. [T8l 
However with larger opacities, optically thick cooling af¬ 
fects larger disk radii, slowing the cooling. Remarkably, 
the adopted MMSN value of opacity hits the a sweet spot 
where optically thin cooling is fast enough, yet slower op¬ 
tically thick can be restricted to inner disk radii. 

At smaller disk radii, Fig. [T8l shows the hallmark of dif¬ 
fusive cooling, that larger wavenumbers can cool faster, 
but not below the floor set by optically thin cooling. 
Optically thick cooling times also rise sharply toward 
smaller radii (as P oc densities and 

short orbital times. This effect suppresses VSI growth at 
small radii, but with a strong wavenumber dependence. 

Fig. [T9| highlights this wavenumber dependence by 
plotting the wavenumbers for which P = Pent- Above 
the solid curves (i.e. for larger wavenumbers), P < Pent 
at all disk heights, implying linear VSI growth. Below the 
dashed curves, VSI growth is strongly suppressed since 
P > Pent for all |z| < 3i7. Between the solid and dashed 
curves some growth is possible, but only fairly close to 
the solid curve (according to >17.31) . Thus linear growth 
of the VSI near 1 AU is only possible if kx P/ lO^/iJ. As 
argued elsewhere, such small-scale disturbances may not 
drive significant turbulence or transport. 

We thus doubt that the VSI is significant at 1 AU in 
MMSN-like PPDs, for any opacity. At higher opacities, 
the required wavenumbers become shorter and even more 
problematic, as shown by the red curves in Fig. [T^l At 
lower opacities, even optically thin cooling (the fastest 
possible) is too slow, as discussed above. 

7.3. VSI growth in the MMSN 

We now consider the growth of the VSI in a MMSN 
disk model. This numerical calculation is similar to ([Hi 
but with different disk parameters and with the self- 
consistent cooling times of Eq. [59] 

We consider the growth timescales of the fundamental 
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Figure 18. Dimensionless thermal relaxation timescales /3, evaluated at the midplane (left) and at z = (right) in the fiducial PPD. 
Eq. 1591 is plotted for three values of the perturbation radial wavenumber: k = 1 (dotted), fc = 10 (dashed) and k = 100 (dashed-dot), for 
three values of the opacity relative to the MMSN: = 0.1 (top), = 1 (middle) and = 10 (bottom). The solid line is the critical 

thermal relaxation timescale dcrit- 



r/AU 

Figure 19. The minimum perturbation wavenumber k in the 
MMSN such that the associated dimensionless thermal relaxation 
time /3 at z = 0 (solid) and z = 3/f (dashed) is less than the criti¬ 
cal timescale / 3crit • Asterisks correspond to the inner cut-off radii 
shown in Fig. 1201 

mode for a range of wavenumbers and disk radii. We 
focus on the fundamental, i.e. lowest order vertical, mode 
because it is the fastest growing mode except for some 
surface modes at high wavenumbers. We neglect these 
surface modes for reasons discussed in ^ and 118.11 


Fig. [201 shows that in the MMSN, the VSI is active 
from CAu ~ 5 to tau ^ 50 with growth timescales ~ 
30—40 orbits. A small radius cut-off exists, inside of 
which growth is strongly suppressed. The cut-off occurs 
at smaller radii for larger wavenumbers, as expected from 
Fig. [TOl Growth at smaller disk radii is possible for yet 
larger wavenumbers, but we remain skeptical about the 
non-linear significance of such short lengthscales. 

The growth times rise gradually as r increases towards 
lOOAU in Fig. [201 This trend is expected as the outer 
radius cutoff at ^ 150AU (from Fig. |TH1) is approached. 
Our numerical results suggest the VSI is efficient at radii 
of few tens of AU i n the MMSN, consistent with esti¬ 
mates made by IN13I using the same disk model. 

8 . CAVEATS AND NEGLECTED EFFECTS 
8 .1. Viscosity 

Our neglect of viscosity is valid in a laminar disk, be¬ 
cause molecular viscosity is so small. However, turbu¬ 
lence could act as an effective viscosity which preferen¬ 
tially damps small scale modes. Since the goal of VSI is 
to drive turbulent transport, the most relevant modes for 
sustaining the VSI should be able to overcome turbulent 
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Figure 20. VSI growth times (tgrow = I/o') in orbital units 
(t/jrb = 27r/r2K) for the reference MMSN disk model with self- 
consistent dust cooling. 

damping. 

To roughly estimate the wavelengths which are prone 
to damping, we consider the stan dard prescription for 
the k inematic viscosity v = aCgH (iShakura fc SunvaevI 
[Ml), in terms of the dimensionless a parameter. The 
viscous timescale for a perturbation lengthscale I ^ 
l/kx = H/k is tvisc = l/aPfl. For significant growth, 
tvisc should be longer than the characteristic VSI growth 
timescale l/{h\q\il), i.e. < \ a\h/a . With h ^ 0.05, 
q = —1 and a from simulations ( N13D . we estimate that 
growth requires fc < 20 for a ~ 10“^ or fc < 70 for 
a - 10 - 5 . 

This argument justifies our focus on moderate 
wavenumbers, k = 0(10). Future work should consider 
the viscous effects in more detail, to better understand 
how the VSI operates in nature and in simulations. 


8 .2. Convective overstability 

The VSI bears s ome similarity to the ‘convective over - 
stability’ ICONO. iKlahr fc Hubbard! 12014 lLvrall2014[) . 
Both instabilities rely on thermal relaxation to overcome 
the stability of disks to adiabatic perturbations, i.e. the 
Solberg-Hoiland criteria (see ^3.41) . 

The key differences are that the CONO neglects verti¬ 
cal buoyancy and requires an imaginary horizontal buoy¬ 
ancy. 


= -■ 


1 dPdS 
pCp dr dr 


< 0 . 


(62) 


Since vertical buoyancy is an important stabilizing influ¬ 
ence for the VSI, it should be considered in future studies 
of CONO. 

The sign of depends on disk parameters. To 
demonstrate that the parameters needed for iV^ < 0 are 
somewhat extreme, or at least non-standard, we consider 
the midplane of vertically isothermal disks with E oc r" 
so that n = p-\- q/2 + 3/2. Using Eqs. I22al and (621 with 
7 = 1.4, V/ < 0 requires n > 3(1/2 -|- q). For stan¬ 
dard disk temperature laws, this requirement implies a 
flat or rising surface density profile, e.g. n = 3/14 for 
q = —3/7 or n = 0 for q = —1/2. Since standard E pro¬ 
files decline with radius, CONO is most like to operate 
at special locations, like the outside edges of disk gaps or 
holes and/or shadowed regions where q < —1/2. 

Perhaps due to these physical differences, CONO op¬ 
erates in different regimes of parameter space than the 


VSI. The CONO grows best for longer wavelengths fc < 1 
and longer cooling times /3 ~ 1. 

Future work should aim to understand these related 
instabilities in the same framework, thereby explaining 
the key differences. 

8.3. Radiative transfer 

Our relatively simple treatment of cooling with an ide¬ 
alized dust opacity could certainly be generalized in fu¬ 
ture works. In hotter disks, gas phase opacities must be 
considered. In cold disks, there are many choices for the 
dust opacity, which varies with grain sizes and composi¬ 
tions. Changing dust properties would alter the viability 
of the VSI for better or worse. The radiative transfer 
itself could be calculated with higher levels of sophisti¬ 
cation, as is already being don e in numerical simulations 
of the VSI (IStoll fc KleviMl . 

9. SUMMARY AND DISCUSSION 

In this paper we study the vertical shear instability 
(VSI) with a focus the role of radiative cooling. In turn 
we assess the viability of the VSI as an angular mo¬ 
mentum transport mechanism in protoplanetary disks 
(PPDs). 

Our linear, axisymmetric analysis of the VSI considers 
(uniquely to our knowledge) finite cooling times in a ver¬ 
tically global model. In order for vertical shear to drive 
the VSI, short cooling times are needed to weaken the 
stabilizing effects of vertical buoyancy. Our main analyt¬ 
ical finding, which we confirm numerically, is the critical 
cooling timescale above which VSI growth is suppressed, 
Eas.[5land[46l 

Our main focus is irradiated, vertically isothermal 
disks which have strong vertical buoyancy. The criti¬ 
cal cooling time is thus short, shorter than the orbital 
time by a factor of the disk aspect ratio. This finding 
is consistent with, and help s exp lain, the results of re¬ 
cent numerical simulations (lN13r ). We briefly consider 
vertically non-isothermal disks in 96.41 

In applying our results to PPDs, we pay particular 
attention to the transition from optically thick radia¬ 
tive diffusion to optically thin Newtonian cooling. The 
largest obstacle to VSI occurs in high density inner disk, 
< 1 - 5AU, where radiative diffusion times are slow. 
Shorter wavelength disturbances cool faster, so long as 
they remain optically thick. In the inner disk, how¬ 
ever, our VSI cooling criterion requires wavelengths that 
are too short to drive significant transport. The best 
hope for the VSI in inner disk is a low surface density, 
which speeds radiative diffusion by lengthening the pho¬ 
ton mean free path. This option naturally begs the ques¬ 
tion of what accretion mechanism lowered the surface 
density in the first place. 

In the outer disk, the VSI tends to cool in the opti¬ 
cally thin limit. The main issue is whether the opacity is 
high enough for optically thin cooling to be sufficiently 
fast. Our standard opacity assumes a Solar abundances 
of small dust. In this case, the VSI can operate from 
^ 5 - lOOAU. A factor of 10 reduction in the opacity, 
for instance by l ocking small dust into pl anetesimals and 
planetary cores (|Youdin fc Kenvonll2013D . makes cooling 
times too slow for VSI growth. In this case cooling is too 
slow not just in the outer disk, but into < 1 AU. 
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An enhanced opacity makes optically thin cooling 
faster and radiative diffusion slower. This shift favors 
VSI growth in the outer disk at the expense of the inner 
disk, pushing the inner limit of VSI growth further out. 
The standard choice — corresponding to Solar abun¬ 
dances in a standard MMSN disk — allows the VSI to 
grow over the widest range of relevant disk radii. 

Our detailed study of the spectrum of VSI modes con¬ 


firms that some artificial ‘surf ace modes’ a re triggered by 
imposed vertical boundaries (|N13l : IBLISD . Domain size 
convergence studies are thus essential. Fortunately, our 
results show that longer cooling times stifle the growth 
of surface modes. Thus, at least in some cases, more re¬ 
alistic radiative transfer also produces more reliable dy¬ 
namics. 

The VSI deserves further study as a viable mechanism 
to drive at least low levels of accretion in cold disks. 


APPENDIX 


A. DERIVATION OF THE APPROXIMATE EQUATIONS 

Here we detail the derivation of Eq. EHl used in the analytical discussion of ^ Starting with Eqs. I27al — I27el we 
set r = I and C = 0 for the vertically isothermal limit and the fully-radially-local approximation, respectively. We 
eliminate the horizontal velocity perturbations {Svx-, Svy) to obtain 


where 


dW ^ dlnp 
dz 


= ^ + (W-Q), 

= iv 

c* D \ D oz 

lU t 

v'^W — 'yv^Q H- {W — Q) = ivc‘, 

ir. 


(? In p \ d5vz 

—— dv^ - m——, 
oz ) dz 


.dlnp 

dz 


(7- l)6v^, 


D = . 


(Ala) 

(Alb) 

(Ale) 

(A2) 


Reduction to a single ODE requires dzD = dzH^- At this point we could apply the low frequency and Keplerian 
approximations to set D ^ ^ then D is vertically constant, and we can obtain Eq. [^Hlmore directly. However, 

to demonstrate that the order of approximation is irrelevant, we will retain dzK^ initially. Using Eq. [T9l—I7T1 


dz dz dr dz dz ) dz 


(A3) 


The function F increases monot onic ally from F = — I atz = 0 to F —)► 2 as |z| 
We eliminate W, Q from Eos. lAIl using Eqs. ISTl and lASl to obtain 


0 =- 


d'^Svz 


dz^ 


1 + 


ik. 


'xC%q 


qclF 


Dr 


+ < P ^ 


kl 

D 


{klcl + xD) Dr^ 
ikxC^q^ In p 

Dr J dz^ 


X' 


dlnp dSvz 
dz dz 

_ c| / dlnp V 
D\ dz j 


ki- 


\kxq 


00, so |F| = 0(1). 


(i-x) + 


X 


qclF 


{klcl + xD) r2 


6vz, 

(A4) 


where we have replaced d/dz by d/dz for the fully-radially-local treatment. Now we make the low frequency and 
Keplerian approximations, setting D —>■ to give 


5v'l 


I +ihqk — 


(k'^ + x) 


qh^F 


In p'Sv'z; + i (x + Inp" — In — ihqk^ 


1-X + 


X 


(k^ + x) 


qh^F 


SVz 


=-v'^ (jt'^ + x'j Svz, (A5) 

in terms of dimensionless variables of Eq. [28l Retaining terms to first order in the disk aspect ratio, d <C I, gives 


0 = Sv'^ + ^I + ihqk^ Inp'Sv'z. + (jt^ + x) + (x + ^hqk^ Inp" — Inp'^ — ihqk^ (I — y) 


Svz 


(A6) 


Approximating the density field by Eq. II7bl then gives Eq. [521 
We can establish a correspondence between our Eq. and Eq. 41 in iLubow fc Pringlel (|I993fl . which describes 
adiabatic axisymmetric waves in a vertically isothermal disk without vertical shear. Accounting for the required change 
of variables, the correspondence is exact after setting q = 0 (no vertical shear) and X = 1/7 (adiabatic flow) in our 
Eq. and applying the approximations in our 114.21 to iLubow fc Pringli . 
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In the fully-radially-local approximation, background radial gradients are ignored except where it appears implicitly 
in the expression for the vertical shear rate (via Eq. HH). This is done by neglecting the terms proportional to C 
in Eq. I27all27bl and l27el i.e. setting C = 0- (Nominally C = 1) which is used in our numerical study.) 

For a power-law disk, the neglected radial gradients are 0(r“^), and they appear in comparison with terms of 0{kx)- 
The neglected terms therefore have a relative magnitude oi 0{h/k), which is small for thin disks (h <C 1) and/or small 
radial wavelengths (fc ^ 1). We show in the following sections that the fully-radially-local approximation only becomes 
invalid in the adiabatic limit, which is not the relevant regime for the VSI. 

We co mment that this approximation is equivalent to that adopted in the vertically global shearing box formalism 
(VGSB, IMP 14) ■ which is an extension of the standard shearing box (|Goldreich fc Lvnden-Belllll965f) to background 
shear flows that are height dependent. 


B.l. Spurious growth of adiabatic perturbations when C = 0 
A limitation of the reduced model described in Appendix ^4.21 and used in (JSJ is that it cannot be employed for 
adiabatic flow when there is vertical shear, even if h/fc <C 1. We explain this by setting /3 —>■ oo and hence y = I /7 in 
Eq. IA 6 l to give 


0 = 5v'f -I- ^1 -I- \hql^ {\np'5vz)' + + k 


7-1 

7 


In p" + P 




Svz. 


(Bl) 


We multiply Eq. IBII by p5v* and integrate vertically, assume boundary terms vanish when integrating by parts, to 
obtain 



p\5vz\'^dz 


i^) rPM\"dz+- rhpSvzYl^dz 

\ 7 J Jzi 7 Jz^ P 

-I-^ P ^1 —^ p\np'^\5vz\'^dz+ \hqk j Itl p'{p5v*)'Svzdz. 


(B2) 


In the presence of vertical shear g ^ 0, Eq. IB2I shows that iP' is complex for real k, which indicates instability for 
any value of 7 > 1. This stability condition is inconsistent with the second Solberg-Hoiland criterion (Eq. [25]) . which 
states that instability requires the disk to be close to neutral stratification (i.e. I 7 — 1 | 1 for a vertically isothermal 

disk). This spurious growth in the C = 0 model arises because we have retained the global radial temperature gradient 
to obtain vertical shear, but have ignored it elsewhere in the linear problem (as well as the background radial density 
gradient). Nevertheless, we demonstrate below that this inconsistency is unimportant for the VSI, which occurs for 
,5 <C I, not for adiabatic perturbations. 


B.2. Effect of global radial gradients 

In Fig. [211 we plot the effect of the global radial gradient terms proportional to ( in Eq. I27ap - I?7el by calculating 
the fundamental VSI growth rates using three approaches. We compute growth rates from the dispersion relation Eq. 
|44l which assumes C = 0; from Eq. I27al -- (27^ with C = 0; and from Eq. I27al — [27^ with C = 1- 

All three methods give similar behavior, and growth rates are in close agreement for fd ^ 1. Differences arise for 
/3 > 1 , and as /3 —>■ 00 the fully-radially-local approximation, where C = 0, gives a (spurious) growth rate as expected 
from the discussion above. Inclusion of the global radial gradient terms results in the expected behavior (cr —>■ 0 as 
/3 —>■ 00 ). Despite this caveat. Fig. shows that provided we consider /3 < 0(1), then setting C = 0 does not affect 
the VSI significantly. 


C. LINEAR PROBLEM IN THE ISOTHERMAL LIMIT 


Here we summarize selected results for isothermal perturbations {/d = 0) in vertically isothermal disks (T 
the fully-radially local approximation = 0). In this case Eq. lAlcI becomes 


Q = W 


1 ) in 

(Cl) 


(i.e. 5P = c^S p). For isoth ermal perturbations it is simpler to work with an equation for W by substituting Q = W 
into Eq. lAlbl and using Eq. lAlal to eliminate Svz- In this case, we shall not yet make the low frequency approximation, 
but first make the Keplerian approximation. We obtain, in terms of dimensionless variables. 


Q = W" + 


In p' 


\hqkf{z) 
1 — 1)2 


W + v^ 



w, 


(C2) 


where for discussion purposes we have defined f{z) such that 


dVL^ 

dz 


h'^qf{z)n'^. 


(C3) 
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Figure 21. Growth rate of the fundamental VS I mode as a function of the thermal relaxation time /3. The ‘Analytical ((^ = 0)’ curv e 
is ca lculated from the dispersion relation Eq. 1441 The ‘Numerical ((J = 0)’ curve is obtained from the numerical solution to Eq. I27al — 
EZH with C = 0, which neglects explicit dependencies on the radial disk structure. The ‘Numerical {C, = 1)’ curve is obtained from Eq. 
[27a]-[27e| with (^ = 1, which accounts for the radial disk structure self-consistently. The disk parameters are ( 7 , T) = (1.4,1.011) and 
(p, O', h) = (—1.5, —1, 0.05). The perturbation wavenumber Is k = 30. 


By comparing Eq. IC3l with Eq. we see that f{z) = z(l + h?‘z^') ^. More generally, though, / may be regarded as 

a representation of the vertical shear profile. Physically, we expect there is a maximum value of \dU,/dz\^ the existence 
of which should limit the growth rate of the VSI, as remarked in ^13.2.21 We explicitly demonstrate this below. 


C.l. Maximum growth rate in the low-frequency limit 

Here we consider the low-frequency limit |z)| <C 1 and show that the growth rate is limited by the maximnm vertical 
shear rate in the domain. We approximate Eq. IC2I as 


0 = W" + 


In p — ihqkf{z) W' -|- (l -f P) W. 


(C4) 


We multiply Eq. IC4I bv pW* and integrate vertically from z = zi to z = Z 2 - We neglect boundary terms when 
integrating by parts, by assuming W or W' vanishes at the boundaries, or that the boundaries are sufficiently far away 
so that the boundary terms are negligible because of the decaying background density with increasing height. Then, 


i? (l -f P) / ' p\W\^dz = f ' p\W'\‘^dz + ihqk f ' pf{z)W*W'dz. (C5) 

It follows that for instability (Imd >0), it is necessary to have q^Q or more generally d^l/dz ^ 0, i.e. there must be 

vertical shear. _ 

The real and imaginary parts of Eq. IC5I are 


(l -k p) y ' p\w\^dz - J ' p\W'\‘^dz = Re 

rZ2 

ihkq / pf{z)W*W'dz 
J Zi 

pZ2 

2u}d (l -k P) J p\W\'^dz = Im 

PZ2 

ihkq / pf{z)W*W'dz 
J Zi 


(C6) 
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where we recall Cj = Ret) and a = Imt). Adding the square of these equations give 

Z2 22 1 ^ 22 

Ivl"^ (l + kA f " p\W\‘^dz- f %\W'fdz +4^2 (l + pW / p\W'\"dz = 

J Z^ J Z^ 'J Z\ 'J Z\ 


It is clear that 


4ct2(i + /c2) / p\W\^dz p\W'\ dz< 

f Z\ J Zi 


hkq / pf[z)W*W'dz 
J Zi 


ihkq / pf{z)W*W'dz 

J Zi 

(C7) 

(C8) 


On the left hand side of this inequality, we apply the Cauchy-Schwarz inequality to obtain 


^ ri2 \ ^22 ^22 

/ p|IV||IVV^ </ p\Wfdz p\W'fdz. 

I 2i } J Z^ J Z^ 


On the right hand side of Eq. IC8I we have 


pf{z)W*W'dz 


< f " p\fiz)W*W'\dz < max (I/I) f %\W\\W'\dz, 
J Zi J Zi 


where max(|/|) is the maximum value of |/| in £ € [zi, Z 2 ]. Inserting these inequalities into Eq. IC8I gives 

|(t| < 


-^lM=max(|/|) < ^max(|/|). 
2 V 1 + A :2 2 


(09) 


(CIO) 


(Oil) 


It follows that the maximum possible growth rate of unstable modes, satisfying the above boundary conditions, is 
limited by the maximum vertical shear rate in the domain considered. 


kr < max 


dfl 

dz 


( 012 ) 


as expected on physical grounds. Note that if the thin-disk approximation is used in an infinite domain, then the 
maximum growth rate is unbounded (since in that case max|/| —>■ 00 ). However, large growth rates invalidate the 
low-frequency approximation and the above analysis breaks down. 

In practice, one might consider a vertical domain of a few scale heights in a thin disk with |g| = 0(1). Then / ~ z, 
so that max|/| = 0(1), implying a maximum growth rate O(hHK), consistent with numerical results. 


C.2. Explicit solutions in the thin-disk limit 

Here we assume a thin disk [h <C 1) so that Inp ~ —z^/2 and /(z) ~ z. However, we do not assume low frequency 
from the onset. Then Eq. 1021 becomes 


W" 



zW' + 



W = 0. 


(013) 


We remark t hat ta king the low frequency limit and considering ^ 1, Eq. l013l becomes equivalent to Eq. 39 in mi 
or Eq. 28 in lBLlSl although we have ta ken a different route. 

We seek power series solutions to Eq. 10131 

00 

W(z)=^a;z'. (014) 

1=0 


Then the coefficients ai are given by the recurrence relation 


{I -|- 2 )(Z -|- l)a ;+2 + 




(015) 


We demand the series to terminate ai I = i.e. a polynomial of order L, so that the vertical kinetic energy density 
remain bounded as |z| —^ 00 . Then the eigenfrequency v is given via 

— (^L 1 L (^1 ihqit^ = 0. (C16) 

Note that we have applied a regularity condition at infinity, since the vertically isothermal disk has no surface. If 
vertical boundaries are imposed at finite height, as done in numerical calculations, then the above solution needs to be 
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modified to match the desired boundary conditions. This enables the ‘surface modes’ seen in numerical calculations 

(IblUi . 


C.2.1. Stability without vertical shear 
In the absence of vertical shear q = 0, Eq. IClOl can be written as 

= (L-v^) (l-v^) , (C17) 

which is just the di spersion relation for axis y mmetric isothermal waves in a yertica ll y isothermal disk fe.g . 
Okazaki_^£_^ Il987t iTakeuchi fc Mivamal 119981 iTanaka et al.l 120021 IZhang fc Tail 120061 IQgilvie fc Latted 120131 
Barker fc Qgilviel 12014 iBLlSf) . In this case the solutions are Hermite poly nomia ls, W oc HeL(z). The eigenfre - 
quency b = d) is real and the disk is stable. The low frequency branch of Eq. IC17I are inertial waves (lBalbusll2003f) . 
For |w| <C 1 and L > 1 the dispersion relation i s = L, or d) oc k~^ for fixed L. This inverse relation has been 

qualitatively observed in numerical simulations of iStoll &: KlevI (1201411 . Note that L = M + 1 , where M > 0 is the mode 
number used for analytical discussion ([JS]) based on the reduced equation for Svz-, rather than for W as considered 
here. 


C.2.2. Instability with vertical shear 

The VSI corresponds to unstable inertial waves. This is readily obtained for large by balancing the last two terms 
of Eq. IC16I to give the low frequency branch. Then 


v^c^L 


1 + iqhk 
L + 1 + 


(CIS) 


which is equivalent to Eq. 34 of lBLlSl in the limit k^ ^ L, 1. This signifies instability for L > 1 since we can choose 
the sign of the square root to make Imb > 0. These are the low frequency unstable modes seen in Fig. |T]and Fig. [71 
for which cr oc |a;|. 

D. ANALYTIC DISPERSION RELATION WITH THERMAL RELAXATION 

D.l. Coefficients 

The coefficients of the dispersion relation Eq. [44] are: 
co = M(M+1)A2, 

ci=i/ 3 |(l- 7 ) l + k^ (l + 2Mf-A\hqkM{M+ 1) + 1)^ , 


C2 = + 1^ A +/ 3 ^ |(1 — 7) 1 + 7fc^(l + 2M)^ — 4ift,gfc7M(M + 1) — 7^A^M(M + 1)| 


C3 = , 


/3|/igfc + 7 i + /igfc^l + 2p^ 


— 3i — 2ifc' 




C 4 = /3^ ^1 + 7 ^^^ 7^1“ ihql^ —2 — ^1 + 
C 5 = 2i/3 ^1 + ^1 + 7p^ , 

C6 = /32(l + 7p)". 


(Dla) 

(Dlb) 

(Die) 

(Did) 

(Die) 

(Dlf) 

(Dig) 


D.2. Finding marginal stability 

To investigate marginal stability we set b = 0 so the frequency, b = d)c, is real; and [3 = the cooling time for 
marginal stability. We take the short radial wavelength limit, ^ 1, of the coefficients in Eq. IDll We consider 
M < 0(1) and assume Pc ^ k- The real and imaginary parts of Eq. 03] then give relations for Cjc and Pc as 

Q=M{M Fl){l-h’^q^k^)F3pchqM{M Fl){Cjck) + [l + 7^^ (1 - 7) (1 + 2M)2] {cbckf (D2) 

+ 2hq'ypc (Cck^ - (ffick^ + Pll'^djl (ffick^ , 

0 =2hqM{M + l)k + pc{l — 7)(1 + 2M)^b ^d)cfc^ + hqk (^Cckj — 2PcUJc (^Cckj — hq^'^Plwc (^Cckj 


(D3) 
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We recall that in the low frequency and thin-disk approximations, h, |wc| ^ 1- We note that |q| = 0(1) and 7 > 1 but 
is 0(1). We assume \cSck\ = 0(1 ), sin c e for inertial waves \uik\ ~ \/l + M. Finally, we further assume /3c -C 1, to be 
justified a posteriori, to give Eq. I45al — l4dbl 


D.3. Maximum critical cooling time 

Here we show that for sufficiently thin disks the longest critical cooling time is that for the M = 0 or fundamental 
mode. This allows us to focus on the fundamental mode to obtain an overall cooling requirement for the VSI. 

Consider the simplified dispersion relations for marginal stability, Eq. I45al — 14^ We write X = uik, 9 = (hqk)^ and 
treat M as a continuous variable. We find from Eq. I45al that 


2X 


dX X^{l-e){2M + 1) 
Zm “ X^ + 2M{M + 9)' 


and from Eq. I45bl that 


2X 


" = [.Y^ + 2M(Af+l)l(^h« 


\ dM 2M + 1 


+ 


d\nX\ 
dM J 


- 2{2M + 1). 


We eliminate dX/dM, making use of Eq. I45al in the process, to obtain 

^dln^c _ il-9)M{M + 1) [2(2M +1)2+ 12X2] (3 2 


(2M + 1)- 


dM 


Hence, 


2 [X2 + 2M{M + 1)(1 - 0)] [X2 + 2M{M + 1)] 


dPc 

< 0 

dM ’ 


(D4) 


(D5) 


(D6) 


(D7) 


for all M > 0 if 0 < 1, which imply max(/3c) occurs at M = 0. This conclusion may also be reached by explicitly 
solving Eq. I45al — l4Fbl with 0 as a small parameter. For fixed k the condition 0 < 1 can be satisfied for sufficiently 
small h. All such modes are stabilized if ,5 > Pc{M = 0) = Pent- 
This result highlights the importance of the fundamental mode — it is the most difficult mode to stabilize with finite 
cooling. Furthermore, for M = 0 we may obtain the expression for /3crit from the dispersion relations Eq. ID21 — |D3] 
without assuming it is much less than unity at the outset or place restrictions on 0. 


REFERENCES 


Armitage, P. J. 2011, ARA&A, 49, 195 

Bai, X.-N. 2015, ApJ, 798, 84 

Balbus, S. A. 2003, ARA&A, 41, 555 

Balbus, S. A., &; Hawley, J. F. 1991, ApJ, 376, 214 

Balbus, S. A., Hawley, J. F., & Stone, J. M. 1996, ApJ, 467, 76 

Barker, A. J., & Latter, H. N. 2015, ArXiv e-prints, arXiv: 1503.06953 

Barker, A. J., & Ogilvie, G. L 2014, MNRAS, 445, 2637 

Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 

Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28 

Blaes, O. M., &; Balbus, S. A. 1994, ApJ, 421, 163 

Chandrasekhar, S. 1961, Hydrodynamic and hydromagnetic stability 

Chiang, E., Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493 
Chiang, E. L, & Goldreich, P. 1997, ApJ, 490, 368 
Fricke, K. 1968, ZAp, 68, 317 
Gammie, C. F. 1996, ApJ, 457, 355 

Goldreich, P., &; Lynden-Bell, D. 1965, MNRAS, 130, 125 
Goldreich, P., & Schubert, G. 1967, ApJ, 150, 571 

Gressel, O., Turner, N. J., Nelson, R. P., Sz McNally, C. P. 2015, ArXiv e-prints, arXiv:1501.05431 
Klahr, H., Sz Hubbard, A. 2014, ApJ, 788, 21 

Lee, A. T., Chiang, E., Asay-Davis, X., Sz Barranco, J. 2010, ApJ, 718, 1367 

Lesniak, M. V., Sz Desch, S. J. 2011, ApJ, 740, 118 

Lesur, G., Kunz, M. W., Sz Fromang, S. 2014, A&A, 566, A56 

Lubow, S. H., Sz Pringle, J. E. 1993, ApJ, 409, 360 

Lyra, W. 2014, ApJ, 789, 77 

McNally, C. P., Sz Pessah, M. E. 2014, ArXiv e-prints, arXiv: 1406.4864 
Nelson, R. P., Gressel, O., Sz Umurhan, O. M. 2013, MNRAS, 435, 2610 
Ogilvie, G. L, Sz Latter, H. N. 2013, MNRAS, 433, 2420 
Okazaki, A. T., Kato, S., Sz Fukue, J. 1987, PASJ, 39, 457 
Salmeron, R., Sz Wardle, M. 2003, MNRAS, 345, 992 
Shakura, N. L, Sz Sunyaev, R. A. 1973, A&A, 24, 337 

Simon, J. B., Bai, X.-N., Stone, J. M., Armitage, P. J., Sz Beckwith, K. 2013, ApJ, 764, 66 

Stoll, M. H. R., Sz Kley, W. 2014, ASzA, 572, A77 

Takeuchi, T., Sz Miyama, S. M. 1998, PASJ, 50, 141 

Tanaka, H., Takeuchi, T., Sz Ward, W. R. 2002, ApJ, 565, 1257 



























22 


M-K. Lin, A. N. Youdin 


Tassoul, J. 1978, Theory of rotating stars 

Townsend, A. A. 1958, Journal of Fluid Mechanics, 4, 361 

Umurhan, O. M., Nelson, R. P., Gressel, O. 2013, in European Physical Journal Web of Conferences, Vol. 46, European Physical 
Journal Web of Conferences, 3003 
Urpin, V. 2003, A&A, 404, 397 

Urpin, V., &: Brandenburg, A. 1998, MNRAS, 294, 399 

Youdin, A. N., & Kenyon, S. J. 2013, From Disks to Planets, ed. T. D. Oswalt, L. M. French, &: P. Kalas, 1 
Youdin, A. N., &; Lithwick, Y. 2007, Icarus, 192, 588 
Youdin, A. N., & Shu, F. H. 2002, ApJ, 580, 494 
Zhang, H., Sz, Lai, D. 2006, MNRAS, 368, 917 



